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Abstract 

We analyze the discontinuous Petrov-Galerkin (DPG) method with optimal test func- 
tions when applied to solve the Reissner-Mindlin model of plate bending. We prove that 
the hybrid variational formulation underlying the DPG method is well-posed (stable) 
with a thickness-dependent constant in a norm encompassing the L2-norms of the bend- 
ing moment, the shear force, the transverse deflection and the rotation vector. We then 
construct a numerical solution scheme based on quadrilateral scalar and vector finite ele- 
ments of degree p. We show that for affine meshes the discretization inherits the stability 
of the continuous formulation provided that the optimal test functions are approximated 
by polynomials of degree p + 3. We prove a theoretical error estimate in terms of the mesh 
size h and polynomial degree p and demonstrate numerical convergence on affine as well 
as non- affine mesh sequences. 

Keywords: plate bending; finite element method; discontinuous Petrov-Galerkin; discrete 
stability; optimal test functions; error estimates 
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1 Introduction 



Finite element methods based on the principle of virtual displacements are the most widely 
used tools for computing the deformations and stresses of elastic bodies under external loads. 
However, in the modelling of thin-walled structures, the basic formulation leads to so-called 
locking, or numerical over-stiffness, unless special techniques (reduced integration, nonconform- 
ing elements) are applied, see [TTj, EU [25]. Another difficulty related to the displacement 
based formulations is the stress recovery. It is well known that the accuracy of the stress field 
derived from the displacement field can be much lower than that of the displacement field. 
Therefore special recovery techniques are often applied to improve the accuracy of stress ap- 
proximations, see [26| 133] . Practical finite element design relies heavily on heuristics, intuition, 
and engineering expertise which make numerical analysis of the formulations difficult, since 
the various physical and geometrical assumptions do not have obvious interpretations in the 
functional analytic setting required for mathematical error analysis. 

Mixed formulations where stresses are declared as independent unknowns are attractive 
because they often avoid the problem of locking by construction and allow direct approximation 
of the quantities of interest. However, in contrast to pure displacement formulations, mixed 
finite element methods do not inherit stability from the continuous formulation, but the stability 
of the discretization must be independently verified for each particular choice of finite element 
spaces as in [21 IH El [131 EEH E321 E0HS2] • The recently introduced discontinuous Petrov-Galerkin 
(DPG) variational framework provides means for automatic computation of test functions that 
guarantee discrete stability for any choice of trial functions, see [T7H2T| I27T429] 135] . 

In this paper we provide an error analysis for the DPG method with optimal test functions 
when applied to the Reissner-Mindlin model of plate bending. We follow the error analysis 
program laid down in [IH1 [23]- The stability analysis utilizes a duality argument based on the 
concept of optimal test space norm and is better suited to multidimensional problems than the 
earlier (see [T8], EH1 ED, E7]) analysis technique based on deriving an explicit expression for the 
generalized energy norm. 

The unknowns in the (mesh-dependent) DPG formulation of the Reissner-Mindlin model are 
the shear force, bending moment, transverse deflection and rotation (field variables) as well as 
their suitable traces defined independently on the mesh skeleton. First, we show that the well- 
posedness and stability of the ideal DPG variational formulation follows from the well-posedness 
of the b ending-moment formulation of the Reissner-Mindlin model which was established in [10], 
see also [U [HI [9] . We introduce then a quadrilateral finite element discretization where the field 
variables are approximated by piecewise polynomial functions of degree p or p + 1 on each 
element and the traces by piecewise polynomials of degree p (resultant tractions) and p + 1 
(displacements) on the mesh skeleton. We prove that on affine meshes the discrete formulation 
is stable in the sense of Babuska and Brezzi provided that the optimal test functions are 
approximated by piecewise polynomials of degree p + 3 on each element. 

The stability estimate is derived using regular (mesh and thickness independent) Sobolev 
norms and the estimate breaks down at the Kirchhoff limit corresponding to vanishing shear 
strains. Our final error bounds are therefore inversely proportional to the slenderness of the 
plate. The analysis indicates that the slenderness dependency arises from the shear stress term. 
This observation is corroborated by the numerical experiments which reveal that the accuracy 
of the shear stress is indeed affected by the value of the thickness while the other quantities are 
rather independent of it. 
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The paper is structured as follows. The derivation of the hybrid ultra-weak variational 
formulation of the Reissner-Mindlin plate bending model is presented in the next Section. The 
wellposedness of the formulation is proved in Section |3j The corresponding finite element 
method is introduced and analyzed in Section [4] and the results of our numerical experiments 
are shown in Section [5] The paper ends with conclusions and suggestions for future work in 
Section |f3 



2 Reissner-Mindlin Plate Bending Model 
2.1 Strong Form 

Let Q be a convex polygonal domain in IR 2 representing the middle surface of a plate. We take 
L = diam(fi) as the length unit and assume that the plate thickness t is small as compared with 
unity, that is the plate is thin. In the Reissner-Mindlin model, the deformation of the plate is 
described in terms of the transverse deflection w and the rotation vector ib, both defined on 
the middle surface fl In the case of linearly elastic, homogeneous, and isotropic material, the 
shear force vector V and the bending moment tensor M are related to the displacements as 
(see for instance V£ 



V = KGtCVw -ib), M = Dt 3 [(l-is)e(tl>) + isti(e('tl>))I], (1) 

where I is the identity tensor and e('0) — ^(^^ + ^^ T ) denotes the symmetric gradient. 
Moreover, 

G=:* d E 



2(l + i/)' 12(l-i/ 2 ) 

are the elastic material parameters written in terms of Young's modulus E and Poisson's 
ratio v while k > is an additional model parameter called the shear correction factor. The 
fundamental balance laws of static equilibrium are 

-V-V = p, -V M-V = 0, (2) 

where p represents a transversal bending load. 
Upon rescaling the static quantities as 

p ■=— >■ Gt 3 p, V m> Gt 3 V, M ■=— >■ Gt 3 M, 

introducing the auxiliary variable uj = |(Vi/> — ^'4 )T )i an d inverting the definition of M in ([T]) 
we arrive at the Reissner-Mindlin system 

k~H 2 V - Vw + ib = 0, C^M - Vib + uj = 0, 

3 

-V-V = p, -V-M-V = 0, V ; 



where 



(TV = 6 I t - ^tr(r)I 



is the two-dimensional "compliance" tensor. 
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2.2 Hybrid Ultra- weak Form 

We use the usual Sobolev spaces H S (X) of scalar- valued functions defined on a domain Icl 2 
and boldface font for the vector- and tensor- valued analogues. As usual, H°(X) = L 2 (X). Ac- 
cordingly, we make use of the space H(div,X) consisting of vector fields in L 2 (X) with diver- 
gence in L 2 (X) and denote by IT(div, X) the corresponding space of tensor- valued functions 
with rows in H(div,X) (the divergence of a tensor is taken row- wise). 

Let {flh} be a non-degenerate family of partitions of Q into convex quadrilaterals, where h 
refers to the maximum element diameter in f^. Integration of the system (|3j) by parts over a 
single element K in fl^ gives 

k-H 2 (V, q) K + (w, V • q) K - (w, q • n) 9K + (ij>, q) K = Vq G if(div, K) 

(C- l M, t) k + (il>, V ■ t) k - (il>, rn) aK + (r J, r) K = VtG H (div, if) 

(V,V^-( 2 ,V.nU = (p,4 V^Gii 1 ^) (4) 

(M, V0) x - (0, Mn) ax - (V, 4>) K = V0 G H l (K) 

(M,sJ)^ = \/seL 2 (tt) 

where n denotes the outward unit normal on OK. The standard L 2 inner product of scalar-, 
vector- or tensor-valued functions over K and dK have been denoted by (•, -)k and (•, -)dK, 
respectively. Moreover, the vorticity has been represented as a single unknown uj = r J, where 

" f 
-1 

and the equilibrium condition M\ 2 = M 2 \ has been imposed weakly using the same notation. 

The next step in developing the DPG formulation is to declare the traces (w, if), Vn, Mn)|g^ 
as indepedent unknowns by rewriting the boundary terms as 

(w, q ■ n) dK ^ (w, q • n) 1/2j9K 

(lf),TIl) dK ^ ($>iTTL)x/ 2 flK 

(z, V ■ n) dK <-> (z, V n )i/ 2j9 K 
(cf>,Mn) dK ^ (0,M„)i/ 2) aK 

where (•, l)x/ 2 ,dK denotes the action of a functional t in YL~ X I 2 acting on scalar- or vector-valued 
functions. 

The boundary conditions for a clamped boundary are ip — 0, w — on dfl and the final 
variational form of the problem is obtained by summing Q over each K in Vt h- The problem 
is to find u = (V, M, w, if), r, w, ip, V n , M n ) G 1A such that 

0(u,v) = £(v) Vv = (q,r,^>, S )GV (5) 

where the functional spaces are defined formally as 



U = L 2 (Q) x L 2 (Q) x L 2 (Q) x L 2 (Q) x L 2 (fi) 

x # 1/2 («9^) x Hl /2 (dn h ) x H-V 2 (dn h ) x H-^idQh) (6) 
V = iZ"(div,fi h ) x Jf (div, n h ) x F 1 ^) x fT 1 ^) x L 2 (Q) 
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and the bilinear and linear forms are given by 

B(u, v) = (V, k-H\ + Vz- 0)„ h + (M, C-V + + s 3)n h 

+ (w, V ■ q)n h + (-0, q + V ■ r)n h + (r J, T)n h - (w, q ■ n}^ - Tn) ffi|i 

- (2, V^)en h - (0, Mn)«i h 
£(v) = (p,z)n h 

Here we have adopted the notation of [23] for elementwise computations of the derivatives on 
the triangulation Qh and its skeleton dQh.: 

(■,-)n h = y~] (•, -)g, (•, •)8n h = (•, -)i/2,dK 

Ken h Ken h 

The broken Sobolev spaces in ([6]) are defined as 

H l {Vt h ) = {ve L 2 (Q) : v\ K E H\K) MK e tt h } 
H(div,n h ) = {qeL 2 (Q) : q\ K e H(div, K) VK e Q h } 

whereas the fractional Sobolev spaces H^ 2 (dQh) and H~ l / 2 (dVt h ) are interpreted as the trace 
spaces of functions in Hq(Q) and H (div, fi) on the skeleton <9i\: 

# 1/2 (<9^) = {v\ dnh : u G #oW 
H- l l\dSl h ) = {77 • n| 9Qh : 77 G Jf (div, fi)} 

The norms in the spaces H^ 2 (dQh) and H~ 1 ^ 2 (dQh) can be defined as 

INI/^fan,) = in /, {IMItfi(n) : 7o(«) = ^} 
\\Vn\\H-V2(dn h ) = M A\\r}\\ H (div,n) ■ Tnfa) = *U 

r)£H (div,S2) 

where 70 and 7 n denote the trace operators satisfying 70(f) = u|an h and 7 n (77) = 77 ■ n\an h f° r 
all f G C x (f2) and 77 G C 1 (f2), respectively. 

3 Well-posedness of the Ultra- Weak Formulation 

We begin with the following formulation of the Babuska-Lax-Milgram theorem and include the 
proof for completeness. 

Theorem 3.1. Assume that 14 and V are two Hilbert spaces and B(u, v) is a bilinear form on 
li x V satisfying 

B(u,v) < C||u|| M ||v|| v VugW,vgV (8) 

sup g | U ' V ^ > «ll v llv Vv G V (9) 
uew Il u llw 

B(u,v) = VvgV u = (10) 
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If C G V' , that is C is a linear functional on V, there exists a unique ueM such that 

B(u,v) = £(v) VvgV 

and 

\\ u \\u < 

a 

Proof. We show that the above assumptions guarantee that also the inf-sup condition 



sup g | U ' V ^ > a\\u\\ u VuGW (11) 
vev H v ||v 

holds. The assertion follows then from the Babuska-Lax-Milgram Theorem, see [21 Theorem 



2.1]. To prove (jn) we define T : U -> V and T* : V ->■ U through 

£(u,v) = (Tu,v) v = (u,TV) w 



It follows from (|8j) and the Riesz Representation Theorem that T and T* are continuous and 
that (19]) is equivalent to 

||T*v|| M > a||v|| v VvgV (12) 

We show next that the range of T* is closed. Namely, if {T*v n } G U is a Cauchy sequence, 
then so is {v„} G V because (|9]) implies that 



vj | v < a\\T*{v m - v n ) I \ u = a\ |T*v m - T*v 



n\\U 



Therefore {v n } converges to some v G V. Because T* is continuous {T*v n } converges to T*v 
which proves that T*(V) = T*(V). 

The condition (fXol) implies now that T* is surjective. If this was not true, there would exist 



a non-zero u£W such £>(u, v) = (u, T*v) = for every v G V. However, this contradicts (10) 



so that we must have T*(V) = ti which together with (12) implies (11): 



B(u,v) (u,T*v) w (u,T*v) w (u,w) M 

sup — : — n — = sup — r: — r; > sup — 77— = a sup — : — r: — = at I |u| \u- 

vev H v ||v vev ll v llv vev a T* v M weu || w ||w 



□ 



3.1 Uniqueness of the Solution 

Lemma 3.1. Let the spaces 14, V and the bilinear form B(u, v) be as defined in Equations ^ 
and 0, respectively. IfuElA satisfies 

<B(u,v) = (13) 

for every v G V, then u = 0. 

Proof. Equation ( [l~3~| ) implies that on every mesh element K we have 

K -1 t 2 (V, q) K + (w, V • q) K - («>, q • n)^ + (ij>, q) K = Vq G ff(div, K) 

(C^M, t) k + (V, V • r)* - rn) ax + (r J, r) x = Vr G H (div, K) 

(V,Vz) x -(z,V;)a K = VzeH\K) (14) 
(M, V<£)* - (0, M„W - (V, 0) x = G H\K) 

(M,sJ) K = \fstL 2 {K) 
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Testing with infinitely differentiable functions which are non-zero only on a compact subset of 
K reveals that 



K~H 2 V - Vw + <0 = 
C^M - Vif> + rJ = 
-V • V = 
-V M-V=0 



(15) 



in every K in the distributional sense. These equations in turn imply that V G U(div, K), 
M G H(dW,K) and w G H\K), if) G H l {K). 
We also have 



w\dK= w\dK, *P\dK= ip\dK and V n \ 9K = V • n\ 9K , M n \ 8K = Mn| S K 



(16) 



This can be seen by integrating each equation in (14) by parts and using the corresponding 
identity in (15) to show that 



(w - w, q • n)i/ 2tdK = Vq G H(dxv, if) 

(V - & Tn) 1/2i9JC = Vr G if(div, if) 

(z,V-n-K) 1/2j ^ = VzeH\K) 

(0,Mn-M n ) 1/2j9K = O Vct>eH\K) 



(17) 



These equations imply that M G i?(div, fl), V G JT(div, fl), and that u> G #o(^)> ^ G ^J(^) 
because w |an= and i/* |sn= 0. 

The extra regularity allows us to set r = M, q = V and z = w, <f> = ?/? in (14). Summing the 
equations together and over every element, we find after integration by parts and simplification 
that 

/t-4 2 (V, V) nh + ( W -us V-n^ + t^ 

(18) 



The second and fourth terms vanish due to ( 17). The last two terms vanish as well. To see this, 
we use ( 16 ) and integrate by parts first locally and then globally (allowed by the regularity of 

(w, V n ) d Q h 



V, w, M, ip) to find that 



(w,V-n) d n h 

(Vw,V) n -(w,V-V)n 
(w, V • n) an 



(19) 



Now the global boundary condition of w G Hq(Q) implies that {w,V n )sn h = 0. A similar 
reasoning and the assumption i/> = Hq(Q) show that (if}, M n )an h = 0. 

Consequently, it follows from ( Jl8| ) that V and M must be zero. To proceed further, we recall 
(see for example [121 Section VI]) that for every r G L 2 (Q), there exists a r r G i?(div, Q) such 
that V ■ r r = and r[ 2 — r r 21 = r. We select r = r r in the second equation of (14) and sum 
over the elements to conclude as in (19) that 



{r,r)n h = (-tp,T r n) 9 n h = (tJ>,T r n) dn = 
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Thus, r = 0. 

Since M and r are already known to vanish, the second equation in (15) implies that ip 
is constant. Since ift G i?g(f2) we find that ip = 0. The first equation in (15) implies then 
similarly that w = 0. Finally ( Jl6| ) shows that also the traces w, ip and V n , M n are zero. Thus, 
all components in u are shown to vanish and the proof is finished. 

□ 

3.2 Existence of the Solution 

In the DPG terminology, the supremum in the condition ^ is called the optimal test space 
norm: 

B(u,v) 



|||v|||v = sup 

ueu \\ u \\u 

In the current application it can be expressed in the form 



\k H\ + Vz-<l)\\l h + \\C V + V0 + sJ|& h + ||V-q||£ h + ||q + V-<r||£ 

\ri2 - T 2 i\\l h + ||[q- n]|||n fc + ||[rn]||| nfe + ||[zn]||| nh + ||[^n]||| nfc , 



h (20) 



where || • ||^ = (•, -)n h and 

lirn nil I Q,m ^ Q ' ll^nlll <nm ^ TJ^h 

||[q ■ n\\\dn h - sup —- , \\[m\\\ 9nh - sup — y- , 

weH^ 2 (8Q h ) W w \\H^(an h ) ^eHl /2 (dn h ) \ m\ H ^(dQ h ) 

\\[za]\\on h = sup -y— , ||[0n]||an h = sup 



It is easy to see that conditions (|9j) and (J8J) of the Babuska-Lax-Milgram Theorem are equivalent 
to the following Lemma. 

Lemma 3.2. There exist positive constants a and C , which are independent of the mesh Vt^, 
such that 

a||v||v < HMHv < C||v||v VvgV. (21) 
Proof. Let v = (q, r, z, <p, s) G V be given and denote by 

(V, M, w, ip, r) G H(div,n) x H(div,fi) x L 2 (ft) x L 2 (f2) x L 2 (Q) 
the solution to the variational problem 

k-H 2 (V, 5V) n + (w,V- 5V) n + (V, <*V)n = (q, 5V) n V<5V G fT(div, fi), 

(C _1 M, 5M) n + (V, V ■ 5M) n + (rJ, 5M) n = (r, 5M) n V<5M G if(div, fi), 

(-V-V.H^^H V5weL 2 (Q), (22) 

(— V ■ M — V, 5if>) n = (0, 5if>) n V<5?/> G L 2 (0), 

(M,5rJ) n = (s,£r) n V5r G L 2 (ft), 

which exists and is unique due to the wellposedness of the bending moment formulation of the 
Reissner-Mindlin model. Namely, the analysis of [10] shows that the bilinear form induced by 



the left hand side of ( 22 ) satisfies the inf-sup condition in a norm encompassing 



*||V||x, a (o), ||V- V\\ L2 (n), ||M|| i2( n), ||V • M + V|| i2(n ), |MUa(n)> \\ip\\ L2 (n), lkl|i 2 (n)- (23) 
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Testing with infinitely smooth functions in the first two equations of (22) reveals then that 
w E H 1 ^), ip G so that the solution of (22) satisfies the estimate 

*HV|U 2(n) + ||V ■ V|U 2(n) + ||M|U a(n) + ||V • M + V|| ia(n) + \\w\\ H i (u) + ||V|| h 1 ^) + I l r l U 2 (n) 

< C (||q|U a (n) + ||r|U a (n) + IMk(c) + ||0IU 2 (fi) + \\ s \\L 2 (n)) 

(24) 

where the constant C > is independent of t, q, r, z, (ft, and s. 



(25) 



The passage from (14) to (15) can be repeated to arrive from (22) to the system 

kTVV - Vw + ip = q 
C _1 M - W + rJ = t 
-V- V = z 
-V-M- V = 4> 

valid on each K in the distributional sense. Now integration by parts yields 

l|q||i 2 (n) + ll T lli 2 (n) + IMI| a (n) + H0lli 2 (n) + ll s llz 2 (n) 

= (k-H 2 V -Vw + il>, q) n + (C^M - + r J, r) n 

- (V • V, z) n - (V ■ M + V, 0) n + (M, s J) n 
= kTV(V, q)^ h + (w, V • q)n h - (w, q • n) a n h + (V>, q)n„ 

+ (C^M, r) Qft + (^, V • r)^ - rn)^ h + (rJ, r) 0h 

+ (V,V*) nfc -<z,V •!!>«,„ 

+ (M, V^)a„ - (0, Mn) fflk - (V, 0) n „ + (M, S J) n , 
Collecting terms and applying Cauchy-Schwarz inequality, we get 



|q|li 2 (n) 



(n) + IkllL(^) + ll ( / > lll 2 (^) + ll s lll 2 (n) 

= (V, «-H 2 q + Vz - cf>) nh + (M, (TV + V0 + S J) n , 
+ (w, V • q)^ + q + V • r)n h + (r J, r)n„ 
- (to, q ■ n)an h - rn) 9slk - V ■ n) af4 - (0, Mn) af!l 



< 1 1 VI 



L 2 



-H\ + Vz- <j>\\ Qh + UMIU^IIC-V + V0 + S J| 



+ lhlU 2 (n)l|V • q||n h + |H|ra(n)||q + V • r||n h + ||r|| ia( n)||ri2 - t 2 i\\l 2 (o) 
+ H[q-n]||anJ|w||ifi(Q) + ll[Tn]|| an J|-i/>||H 1 (fi) 
+ HNI|9nJ|V|| ff(diVj n) + ||[^]||anJ|M|| H(diVi n) 
< 2|||v|||v(HV|| Jf( div,n) + ||M||ir(div,n) + WMln^n) + IMIj/ 1 ^) + IMk(n)) 



By using the estimate (|24|), we obtain 

HI<£lll 2 (n) + IHII 2 (n) 

|v|||v (||q|Ua(n) + \\T\\L 2 (n) + \\z\\l 2 (q) + ||0lk 2 (fi) + IMU 2 (n)) 



|q|li 2( n) + 



\L 2 (Q) 



+ z 



2 

L 2 (Q) 

< cr l \ 



and, consequently, 

l|q|U 2 (n) + IMU 2 (ft) + IMU a (n) + II^IUa(n) + IMk(n) < C^IIMHv- 



(26) 
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The remaining terms constituting the norm ||v||v can be bounded from above by 
directly or by using the triangle inequality: 



Mllv 



< cr^lvlllv 



|V • q||n h < |||v|||v 

|V ■ r||n h < ||V ■ t + q||n h + IMU 2 (n) 
||Vz||n ft < \ \K-H 2 q + V z - (f>\\n h + K-H 2 \\q\\ L2{n) + 
llV^ln^HC- 1 -! 



L 2 (n) 



< cr^iHHv 



(27) 



V0 + sJ||n h + \\C l r\\ L2{ n) + 2\\s\\l 2{n) < Ct ^IMHv 



The first inequality in (21) follows now from (27) and (26) with an a proportional to t. 

The proof of the second inequality is more straightforward. The integral terms || • can 
be bounded from above by ||v||v using the triangle inequality whereas the jump terms can be 
handled by integration by parts and Cauchy-Schwarz inequality: 



W- n \\\dn h 



sup 



(z, q ■ n) 



\ z \ I Hi(Q) zeHl^l) 

Similar arguments can be used to show that 



\ z \\m{tt) 



|q||H(div,n A ) 



IIMHaOh < IM|jT(div,n h ) 
1 1 Milan,, < ||«||flj(n fc ) 

li^iU<ii0il^) 

We leave the details to the reader and conclude our proof. 



□ 



We have shown in Lemmas 3.1 and 3.2 that the conditions of the Babuska-Lax-Milgram 
theorem 3T hold. In other words, we have established 

Theorem 3.2. The ultra-weak variational formulation of the Reissner-Mindlin plate bending 
problem defined by (JsJ) (JtJ) is well-posed. 



Remark 3.1. The proportionality of a to t, in (21), is due to the first term in (24) which 
affects only the shear stress. This observation is ratified in our numerical experiments below. 



4 The Approximate Problem 

In order to discretize (|5|), we choose a finite element trial function space 14 h C U and construct 
a corresponding test function space V r h = T r (lih) C V C V by solving the auxiliary problem 

(T r w h ,v) v =B(w h ,v) VveV r 

for each G Lit. The discontinuous Petrov-Galerkin approximation 6 Uh is defined as the 
solution to the problem 

B(u„,v) = £(v) VvgVI (28) 

The space V r is determined by an appropriate enrichment of the trial function space tih. The 
level of enrichment is specified so that the Fortin's Criterion for the discrete inf-sup condition 
holds: 
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Lemma 4.1. (Fortin's Criterion for DPG) Suppose that for the subspaces V r , Uh, there exists 
a bounded linear projector Ylh '■ V — > V r such that 



B{w h , v - n h v) = \/w h £ U h . (29) 
If 1 1 lift 1 1 < c, then the finite element spaces Uh and V r h satisfy the inf-sup condition 

sup B ^ h r f h) > -\\u h \\ u Vu fc G U h (30) 
vVeVl ||v h || v c 



and £/je DPG approximation is uniquely defined by (28) and is a quasi-optimal approximation 
of u, namely 

Cc 

||u-u h || M < — min ||u-w ft || M (31) 
Proof. See proof of Theorem 2.1 in [23]. □ 



To make Lemma applicable in the present context, we need to construct local projectors 
from ff(div, K) and H X (K) to suitable finite element spaces. In |23j, these projectors were 
constructed for polynomial spaces on simplicial triangulations of Q. We will use the techniques 
of [3] to construct analogous projectors for quadrilateral meshes. We assume the partitions 
to be shape-regular in the usual sense, that is, each angle of each K £ Qh is assumed to be 
bounded away from and n by an absolute, positive constant and the ratio of any two sides 
on K is assumed to be uniformly bounded. 

Let K be a rectangular reference element, and denote by : K — > M 2 the bilinear 
diffeomorphism onto the actual element K = Fk{K)- We define the local bilinear quadrilateral 
finite element space of degree r as 

S r (K) = {ve L 2 (K), v = voF K \ve Q r (K)}, (32) 

where Q r {K) = V rtT {K) denotes the space of polynomials of degree at most r in each variable 
separately on K. We also use the local vector finite element space 

V r (K) = {q : K R 2 | q = (P K q) o ¥~ K \ q £ KT r (K)}, (33) 

where 7V7~ r (K) = V r +\ >r {K) x V r ,r+i(K) is the Raviart-Thomas space and denotes the 
Piola transformation which is defined in terms of the Jacobian matrix 3^ = D^k as 



det Je-(x) 



For the numerical fluxes and traces we need local polynomial spaces defined on the boundary 
dK as 

T r (dK) = {7 £ L 2 (dK), j\ E e V r {E) for all edges E of K}, 
f r (dK) = T r (dK)f]C(dK), 

where V r (dK) stands for polynomials of degree r on E and C(dK) stands for the space of 
continuous functions on dK. 
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The trial space of degree p for the DPG method is defined in terms of the above space^]as 
Uh = {(V, M, to, ij), r, w, xj), V n , M„) G U : 

V\ K e V P (K), M\ K e V P (K), w\ K e S P (K), ^\ K e S P (K), r\ K e S P (K), 
\ 9K e f p+1 (dK), i{>\ dK e f p+1 (dK), V n \ 9K e T p (dK), M n \ 9K E T p (dK) \/K G Q h } 



w . 



In the definition of the enriched test function space V r , we may employ the space (32) 
to approximate those components which belong to H X {K) or L 2 (K) and the space (33) to 
approximate the components in H (div, K). The definition of V r is 

V r = {(q,T,z,(f>,s)eV : q\ K eV r (K), r\ K eV r (K), 

z\ K e S r (K), <f>\ K e S r (K), fi\ K e S r (K) MK g n h }. 

Next we will show, that taking r = p + 3 is sufficient to guarantee the existence of the 
projector needed to guarantee the best approximation property of in Lemma 4^ The proof 
consists of three parts and follows closely the reasoning used in [23] with small modifications. 

Lemma 4.2. Let B(K) be defined as B(K) = {zG S P+2 (K) : z is zero at the vertices of K}. 
Then there exists a projector R° K onto B(K) C H 1 ^) such that 

(R° K z,v) K = (z,v) K VveS p (K) (34) 
{R Q K z,l)dK = {z,i) d K VieT p (dK) (35) 
h~K \ \^°k z \\l 2 (K) + \R%z\ H i^ K ) < C(hx \\z\\ La (K) + \ z \h' l (k)) (36) 
for allz G H^K). 

Proof. To see that R° K is well-defined, we first note that the number of conditions in (34) and 
(35[ is 

dim S P (K) + dimT p (d K) = (p + l) 2 + 4(p + 1) = p 2 + 6p + 5 

and equals the dimension of B(K): 

dimB{K) = (p + 3) 2 - 4 = p 2 + 6p + 5 

Therefore, in order to show that R° K z exists and is unique, it suffices to show that z = implies 
R\z = 0. On each edge e of dK, R\z has the form R° K z\ e = B e u where u G V p (e) and B e is a 
quadratic bubble function defined on e such that < B e < 1. Consequently, (35) implies that 
R° K z\ e = on each edge. This in turn means that R° K z = Bx4>p, where (f> p G Q P (K) and is 
the biquadratic bubble function defined on K such that < Bk < 1 and Bk\ok= 0. Now (34) 
implies that R° K z = 0. The mesh regularity hypothesis and a scaling argument guarantee the 
validity of (36) with a constant C independent of K. □ 

We can now construct a projector into the enriched finite element space such that the 
i7 1 -norm is bounded by an /i-independent number. This is the content of the following Lemma. 

Lemma 4.3. There exists a projector Rk from H 1 ^) into S P+ 2{K) such that 

(R K z,v) K = (z,v) K \fvES p (K) (37) 
(R K zn)dK = (zn)dK VieT p (dK) (38) 
WRRzWmiK) < C|M|ffi(ir) (39) 

for allze H l {K). 



1 A tensor- valued function is included in V P (K) row- wise according to the definition (33 1. 
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Proof. Rrz is defined as Rkz = R%(z — z) + z, where z is the constant function 

LzdK 

which, by a scaling argument and a variant of Friedrichs' inequality, satisfies 



\z — z 



\l 2 {k) < Ch K \z\ H \( K ) 



It follows from the definition of Rk that Rkz — z = R K (z — z) — (z — z) so that (34) and (35) 
imply (37) and (38). 



We have 



\Rkz\\l 2 {k) < \\R K i z - z )\\l 2 {k) + \\z\\l 2 (k) 



□ 



Lemma 4.4. There exists an operator ttk '■ H(div,K) — > V P+ 2{K) such that 

(q-ir K q,r]) K = Vry G S P (K) 
{j, (q - 7r K q) ■ n) dK = V 7 G f p+1 (<9iT) 

ll 7r /<q||//(div,if) < C1|q||jf(div,^) 



(40) 
(41) 
(42) 



/or all q G i?(div, iT). 

Proof. We start by constructing a bounded projector 7r^ : H(div, K) — >■ Qp +2 (i^) for the 
rectangular master element i^. The construction is based on the observation that (40) and (41 ) 
resemble closely the canonical degrees of freedom in the Raviart-Thomas space 1ZT p+ i(K) = 
V p+ 2,p+i{K) x 7^4.1^4.2 (K). Namely, if we denote by T^^dk) the L 2 (9i^)-orthogonal com- 
plement of T p+ i(dK) in T p+1 (dK), and define 

n(k) = { q g nr p+l {k) ■. ( 7 , q • h) dK = v 7 g r^dK)} 

then the operator 7r^ : //"(div,^) — > H{K) is indeed well-defined by the conditions 

feq, »))k = (q, v *7 G V P . P+1 (K) x V p+ i iP (K) 

(TT^q ■ n, 77)^ = (q • n, r?)^ Vr/ G f p+1 (9^) 

This is true because 7r^-q G Tl{k) is a function in lZT p+ i(k) and all of its degrees of freedom 
must vanish when q = 0. 

The corresponding projection for an arbitrary element K = Fk{K) can be defined using 
the Piola transform as tz k = P K o tz ^ o P~ 1 . We have J^r; G V PtP+ i(K) x V p+ i iP (K) whenever 
G Qp{K) so that (40) and (41) follow from the identities 



(q - 7Z K q, r)) K = (Jjf(q - TT^q), f]) k = (q - TT^q, J^t)) 
((q-7r^q) •n,7) sx = ((q - 7T^q) -n,7)aK 



A' 



To prove ( |42| ), we first assume that /i^ = 1 and notice that 7T£ from i?(div, K) to L/2(K), 
Pk from i?(div, K) to H (div, K) and P^ 1 from i3"(div, if) to H (div, K) are bounded opera- 
tors with bounds depending only on the shape of K. Therefore, ttk is bounded from H (div, K) 
to L 2 (K). 
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To extend the ^(-fO-bound to an arbitrary convex quadrilateral K, we follow [3] and 
introduce the dilated element K = D(if) defined by D(x) = h]^x. We have then F^> = DoF A 
so that ir^ = o tz k o Pt 1 and for any q £ iT(div, K), let q = h K q(h K x). Then, 

l|qlL 2 (£) = llqlU 2 (K) 

l|V-q|| L2( ^) = h 2 K \\V ■ q\\ L2{R) = h K \\V -q|U 2 (if) 



so that we have 

IK-K-qlUaCfO = \\ h K 7r K^\\L 2 (K) = 1 K<l\ \l 2 (K) - C \ 1^1 \H(div,K) - C '( 1 1^1 U 2 {K) + h K \ | V-q| | L 

To obtain an /i-independent bound for the norm of the divergence, we use the identities 

(V-q)oF,= ^ 
det Ja-(x) 

V • (TT^q) = Q p+ iV ■ q 
where Q p +i denotes the L 2 (i^)-projector onto Q p+ i(K), to write 

V-(7r^q) ftp+iV-q fl p+1 [detJ^(V-q)oF A -] 



{K)j 



det J A det J A det 3 K 

In other words V • (^i^q) = A A (V • q), where Ax : L 2 (K) — > L 2 (K) is defined by 

fl p+1 [det(JK)(/°F A -)] 

A ^ = 1 ,/ T \ K 

det(J A ) 

for any scalar function /. Now M2j follows because: 

||AK/IU 2W <C||/|| i2(A) V/£L 2 (iT) (43) 



The bound (43) is obvious for elements with unit diameter and can be extended to elements 
with arbitrary diameter with a constant depending only on the shape of K by using the dilation 
x <^-> frj^x. □ 

We can now state our main approximation result: 

Theorem 4.1. Let u = (V, M, w, if), r^, w, -0, V n , M n ) denote the exact solution to the Reissner- 
Mindlin model and Uh = (V/,, Mj, w/,, ^, r/,, 4i Mn,ii) the DPG approximation of 

degree p on an affine mesh with maximum element diameter h. The approximation error 

e = ||V - Vh\\L 2 (n) + ||M - M h \\ L2(n) + \\w - w h \\ L2(n) + - */>JU 2 (n) + \\r - r h \\ L2{n) 

+ \\w- U) h \\ H i/2( dnh) + HV' - V'JIiji/z^) + \\V n - V n ,h\\H-V 2 {dQ. h ) + l|M„ - 
satisfies an a priori estimate 

e < Cr 1 + \\M\\ HP+ 2 {n) + \\w\\ HP+ 2 {n) + \\^\\ HP +2 {n) ) (44) 

where the constant C is independent of h and t but depends on p and Q. 
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Proof. We start by defining a global projection operator 11^ : V — > V p+a piecewise^J 

(II fc v)|jc= (Tr K T,ir K q,R K z,R K ct>,Q K fi) 



where ir^ and Rk are the projectors defined in Lemmas 4.4 and 4.3 and Qk is the L 2 -projector 
onto S P+ 2(K). The projectors satisfy 



(k 1 t 2 r) h + 6 h , q - 7T K q) K 
{C~ x a h + p h 3, (r - 7t k t)) k 



0. 
0. 



(v h , V- (q-7T K q))^ = 0, 
'6 h ,V • (t - ic k t)) k = 0, 

{ilh^( z ~ Rrz)) k = 0, 
(a h ,V(ct>-R K ct>)) K = 0, 



(a h , (s - Q k s)3)k = 



(v h , (q - 77 K q) ■ n) 1/2,9 a- = 0, 
{6 h , (t - 7r K T)n)i /2t dK = 0, 

(Z - R K Z,fln,h)dK = 0, 
(0 — Rk4>, & n)l/2,8K = 0, 

(45) 

for all = (ri h ,a h ,v h ,0 h ,p h ,v h ,e h ,fj n;h ,a- nth ) G The first and third columns follow 

directly from Lemmas 4.4, 4.3 and the definition of Qk- The second column is proved using 
the same Lemmas in conjunction with integration by parts. The first equality in the second 
column holds because 

(v h , V • (q - n K q)) K = (v h , V • (q - ^ K ^))k 

= (Vh, (q - *>q) • fi W - (Vuft, q - Tr K ci)K 
= 

The third equality in the second column holds because 

(rj h , V(z - R K z)) K = (rfh • n, z - R K z)dK - (V • ri h , z - R K z) K 
= (rj h ■ n, z - R K z) dK - (V • f] h , z - R K z) K 
= 

The second and fourth equality can be proven in the same way so that we have established 
the condition (29) of Lemma 4.1 In other words, we have established the best approximation 
property (31 ). 

A more quantitative error estimate is obtained by using results from approximation theory. 
For smooth enough vector and scalar fields V and w, there exist interpolants V G H (div, f2/J 
and w G i/ 1 (fi/ l ) such that 



V\ K eV p {K), w\ K eS p (K) VKeSl h 



and 



IV- VI 



£ 2 (fi) 



\w - vj\\L 2 (n) 



<Ch p+1 \\V\\ HP+1{n) , 
<Ch p+l \\w\\ HP+Hn) . 



Here w is the standard interpolant of w and V denotes the projection of V to the Raviart- 
Thomas space, see [31 122]. 

We can also construct interpolants satisfying w|an h G H^ 2 (dQh) and V-n^^G H^ 1 ^ 2 {dVlh) 
such that 



w 



dK eT p+1 (dK), V-n\ dK eV p (dK) VKen h . 



2 The operator ttk acts on tensors row-wise. 
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Since the traces w and V n associated to the exact solution equal the traces of the corresponding 
field variables, we are allowed to write 

min | \w - v h \\ H i/2 {dnh) < \ \w - w\\m(n) 

min || V n — fj nt h 1 1 #-1/2(0^) < ||V — V||jf(di v ,n) 

Vn,h 

Defining w as the regular interpolant of w with quadrilateral elements of degree p + 1 and V as 
the projection of V into the Arnold-Bofn-Falk space of index p, we obtain the error estimates 

\\w-w\\ H i(n) < Ch p+1 \\w\\ HP +2 {n) , 
||V-V|| ff(div , n) <Ch p+1 (\\V\\ HP+1 + \\V-V\\ 

Identical constructions can be carried out for the remaining solution components ?/>, M, r, 



t/j, and M n . Hence, the estimate (44) is established. □ 



Remark 4.1. Notice that the restriction of the proof to ajfine mesh sequences arises from the 



terms involving r\ h and cr ft in the first column of (45). Namely, when the mapping ¥^ is 
not affine, the use of Piola transform introduces a non- constant factor 1/ det 3k violating the 
orthogonality conditions established in Lemma\4-4\ On an affine mesh, the same terms dictate 



the enrichment degree to be three, since we need to apply Lemma 4-4 a ^ so when r) G V P (K). 
On the other hand, the use of Piola transformation for the shear force and bending moment is 
necessary in general to match the normals in V n and M n with the ones in V • n and Mn. 

Remark 4.2. When bounding the approximation error ofV n and M„ ; use of Raviart- Thomas 
projector would imply loss of one power of h in the convergence rate on a general mesh, see 
|3J/. In the DPG approximation the resultant tractions can be extended as well to the mentioned 
Arnold- Bo ffi-Falk space defined on the reference element as A33J- V [K) = V P+ 2, P (K) xV p , p+ 2(K) 
since the normal components of the elements of this space are also polynomials of degree p on 
the edges. 

5 Numerical Results 

We study the convergence of the DPG method when applied to solve the model problem pro- 
posed in [15] . The problem consists of a fully clamped, homogeneous and isotropic square plate 
loaded by the pressure distribution 

p(x, y) = -— 1 2 J 12y(y - l)(5a; 2 - 5x + l)(2y 2 (y - l) 2 + x(x - l)(5y 2 - 5y + 1)) 
1/(1 — v z ) 

+ I2x(x - l)(5y 2 - 5y + l)(2x 2 (x - l) 2 + y(y - l)(5x 2 - 5x + 1))] 

on the computational domain Q = (0,1) x (0,1). The problem has a closed form analytic 
solution than can be used to address the accuracy of numerical solution schemes. 

We use the values v = 0.3 and k = 5/6 for the Poisson ratio and the shear correction factor, 
respectively. We set p = 1 and compute the DPG solution using uniform and trapezoidal 
iV x iV-meshes, with TV varying as N = 4, 8, 16, 32, 64, see Fig. [TJ 
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(Uniform) 




(Trapezoidal) 



Figure 1: Uniform and trapezoidal mesh sequences. 

The results for the thickness values t = 0.1 and t = 0.001 are summarized in Figs. [2] and 
[3j respectively. In the figures we show the relative errors in the L 2 norm for all quantities of 
interest * 

||V- VfeH HM-M/,11 \\w-w h \\ H^-^JI 

||v|| ' ||M|| ' |H| ' ll^H 

The results show that 

Uniform Trapezoidal 



10° 




Figure 2: Convergence at t = 1/10. Uniform versus trapezoidal mesh sequence. 



1. Optimal quadratic convergence is attained for all quantities on both mesh sequences at 
t = 1/10. 

2. Convergence of the shear stress slows down at t = 1/1000 especially on the trapezoidal 
mesh sequence. However, a relative error of less than 10 percent is attained also at the 
16 x 16 trapezoidal mesh. 
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* — * bending moments 
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10" 
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10" 




shear force 
bending moments 
deflection 
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10"' 
mesh size, h 



Figure 3: Convergence at t = 1/1000. Uniform versus trapezoidal mesh sequence. 




-0.002-0.001 0.001 0.002 
-0.00286 0.00286 



-0.002-0.001 0.001 0.002 
-0.00286 0.00286 



Figure 4: Shear forces at t = 1/1000. 



Finally, we show in Figs. |4j^7] contour plots of all quantities of interest at t = 1/1000 
obtained with DPG by using a fine mesh. The good approximation quality of all quantities 
makes prediction of the values and the locations of maximum stresses straightforward. 



6 Concluding Remarks 

We have analyzed the discontinuous Petrov-Galerkin finite element method in the Reissner- 
Mindlin plate bending problem. The formulation is based on a piecewise polynomial approxi- 
mation using quadrilateral scalar and vector finite elements of degree p for all quantities of in- 
terest (shear stress, bending moment, transverse deflection, rotation). In addition, the resultant 
tractions and the kinematic variables are approximated on the mesh skeleton by polynomials 
of degree p and p + 1, respectively. 

We have showed that the non-standard variational formulation underlying the DPG method 
is well-posed. Based on that result, we have showed that a discretization where the test func- 
tions are approximated in an enriched finite element space of degree p + 3 is stable as well and 
leads to optimal order of convergence in the L 2 norm for all variables. However, the theoretical 
stability estimate breaks down at the limit of zero thickness and therefore the final error bound 
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-0.00028 0.00028 -0.00028 0.00028 



Figure 6: Rotations at t = 1/1000. 
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2e-^^^^^^6e-5 8e-5 
-4.67e-9 8.146e-5 

Figure 7: Transverse deflection at t — 1/1000. 

becomes amplified by the factor t" 1 . Our numerical results indicate that some error amplifi- 
cation indeed occurs for the shear force, but the obtained stress values are relatively accurate 
even on severely distorted meshes. Future work includes formulation of the algorithm for more 
general geometries and an evaluation of the computational cost and robustness in comparison 
with other type of formulations. 
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